This reproducible R Markdown
analysis was created with workflowr (version
1.7.0). The Checks tab describes the reproducibility checks
that were applied when the results were created. The Past
versions tab lists the development history.
Great job! The global environment was empty. Objects defined in the
global environment can affect the analysis in your R Markdown file in
unknown ways. For reproduciblity it’s best to always run the code in an
empty environment.
The command set.seed(20240716) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Using absolute paths to the files within your workflowr project makes it
difficult for you and others to run your code on a different machine.
Change the absolute path(s) below to the suggested relative path(s) to
make your code more reproducible.
Great! You are using Git for version control. Tracking code development
and connecting the code version to the results is critical for
reproducibility.
The results in this page were generated with repository version 18fb11e.
See the Past versions tab to see a history of the changes made
to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Note that any generated files, e.g. HTML, png, CSS, etc., are not
included in this status report because it is ok for generated content to
have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/index.Rmd) and HTML
(docs/index.html) files. If you’ve configured a remote Git
repository (see ?wflow_git_remote), click on the hyperlinks
in the table below to view the files as they were in that past version.
File
Version
Author
Date
Message
Rmd
18fb11e
avm27
2024-07-17
Testing
html
23fe19a
avm27
2024-07-17
Build site.
Rmd
52486bd
avm27
2024-07-17
Testing
html
6284345
avm27
2024-07-17
Build site.
Rmd
98eff28
avm27
2024-07-17
Testing
html
5bc128a
avm27
2024-07-17
Build site.
Rmd
292e3b5
avm27
2024-07-17
Testing
html
12fbe39
avm27
2024-07-17
Build site.
Rmd
479c9dd
avm27
2024-07-17
Testing
html
25a76da
avm27
2024-07-17
Build site.
Rmd
89ef098
avm27
2024-07-17
Testing
html
cb49269
avm27
2024-07-17
Build site.
Rmd
ca44efa
avm27
2024-07-17
Testing
html
21b70e3
avm27
2024-07-17
Build site.
Rmd
ddb0db7
avm27
2024-07-17
Testing
html
a1c090c
avm27
2024-07-17
Build site.
Rmd
adbd778
avm27
2024-07-17
Testing
html
dc04b15
avm27
2024-07-17
Build site.
Rmd
65a97a6
avm27
2024-07-17
Testing
html
e58ce5e
avm27
2024-07-16
Build site.
Rmd
ad8ae96
avm27
2024-07-16
Start workflowr project.
Cocaine craving drives repressive epigenetic
differential gene expression in ventral tegmental area dopamine
neurons
Tate A. Pollock1,2, Alexander V. Margetts1-3, Samara J. Vilca1-2
& Luis M. Tuesta1-3*
1 Department of Psychiatry & Behavioral Sciences 2 Center for
Therapeutic Innovation 3 Sylvester Comprehensive Cancer Center
University of Miami Miller School of Medicine, Miami, FL 33136 *
Corresponding author (ltuesta@miami.edu)
Abstract
Dopamine (DA) signaling plays an essential role in reward valence
attribution and in encoding the reinforcing properties of natural and
artificial rewards. The adaptive responses from midbrain dopamine
neurons to artificial rewards such as drugs of abuse are therefore
important for understanding the development of substance use disorders.
Drug-induced changes in gene expression are one such adaptation that can
determine the activity of dopamine signaling in projection regions of
the brain reward system. One of the major challenges to obtaining this
understanding is the inherent cellular heterogeneity in the brain, where
each neuron population can be defined by a distinct transcriptional
profile. To bridge this gap, we have adapted a virus-based method for
labeling and capture of dopamine nuclei, coupled with nuclear
RNA-sequencing and bioinformatics analyses, to study the transcriptional
adaptations, specifically, of dopamine neurons in the ventral tegmental
area (VTA) during cocaine taking and cocaine craving using a mouse model
of cocaine intravenous self-administration (IVSA). Our results show
significant changes in gene expression across non-drug operant training,
cocaine taking, and cocaine craving, highlighted by an enrichment in
expression of repressive epigenetic modifying enzymes during cocaine
craving. Immunohistochemical validation further revealed an increase of
H3K9me3 deposition in DA neurons during cocaine craving. These results
demonstrate that cocaine-induced transcriptional adaptations in dopamine
neurons vary by phase of self-administration and suggest that such an
approach may be useful in identifying relevant phase-specific molecular
targets to alter the behavioral course of substance use disorders.
Differential Gene Expression Analysis with
DESeq2
0.1 Setup
Building DESeq2 input files from gene count matrices output by
StringTie following removal of samples that are outliers.
Running DeSeq2 analysis with factor levels of FoodTrained vs
Maintenance, Craving vs Maintenance and Craving vs Foodtrained
1 Plot Estimates and
Gather Results
1.1 Gather
Results
These results include overview for each result file and dispersion
estimates based on count values.
1.1.1 Craving vs
FoodTraining
Craving vs FoodTrained CIVSA Results Overview
baseMean
log2FoldChange
lfcSE
stat
pvalue
padj
Mrpl15|Mrpl15
2162.26767
0.3962276
0.3468265
1.1424376
0.2532722
0.9784152
Gm39587|Gm39587
20.30129
-2.1136742
3.3911278
-0.6232954
0.5330904
0.9997970
Lypla1|Lypla1
994.90242
0.0198799
0.8799100
0.0225931
0.9819748
0.9997970
Xkr4|Xkr4
4459.42848
-0.0337912
0.2715766
-0.1244260
0.9009780
0.9997970
Gm39585|Gm39585
5797.38480
0.0861245
0.2414588
0.3566843
0.7213282
0.9997970
Tcea1|Tcea1
2657.79689
-0.4037632
0.3330614
-1.2122786
0.2254058
0.9485680
Rgs20|Rgs20
422.79650
-0.1889366
0.7619761
-0.2479561
0.8041684
0.9997970
Gm16041|Gm16041
21.30106
7.5945314
3.6799019
2.0637864
0.0390380
0.6410945
Atp6v1h|Atp6v1h
10659.10495
-0.4956230
0.2967682
-1.6700675
0.0949060
0.8425404
LOC108167788|LOC108167788
11.33585
-0.8220772
3.4983015
-0.2349933
0.8142140
0.9997970
Npbwr1|Npbwr1
924.87540
-0.3079659
0.6636912
-0.4640199
0.6426335
0.9997970
Oprk1|Oprk1
15661.74234
-0.1442457
0.3231485
-0.4463757
0.6553259
0.9997970
St18|St18
279.69289
-2.0082890
0.9250037
-2.1711146
0.0299225
0.5845690
Rb1cc1|Rb1cc1
12217.86115
-0.0796796
0.2844821
-0.2800866
0.7794110
0.9997970
4732440D04Rik|4732440D04Rik
5369.16278
-0.4918180
0.2427412
-2.0261004
0.0427545
0.6632640
Gm26901|Gm26901
163.46252
-0.2293404
0.9562622
-0.2398300
0.8104620
0.9997970
Pcmtd1|Pcmtd1
4321.05905
0.2034726
0.2645104
0.7692425
0.4417494
0.9997970
Rrs1|Rrs1
625.86245
-0.4850563
0.6605362
-0.7343371
0.4627433
0.9997970
Mybl1|Mybl1
184.22714
1.4451732
1.2066675
1.1976566
0.2310507
0.9561403
Adhfe1|Adhfe1
452.68051
-0.1536903
0.7879557
-0.1950494
0.8453543
0.9997970
out of 23922 with nonzero total read count adjusted p-value < 0.1
LFC > 0 (up) : 250, 1% LFC < 0 (down) : 48, 0.2% outliers [1] :
2272, 9.5% low counts [2] : 1352, 5.7% (mean count < 2) [1] see
‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’
argument of ?results
1.1.2 Maintenance vs
FoodTraining
Craving vs FoodTrained CIVSA Results Overview
baseMean
log2FoldChange
lfcSE
stat
pvalue
padj
Mrpl15|Mrpl15
2162.26767
0.6589546
0.3707556
1.7773286
0.0755142
0.3782100
Gm39587|Gm39587
20.30129
-1.7619611
3.6258386
-0.4859458
0.6270056
0.9345721
Lypla1|Lypla1
994.90242
0.0914200
0.9407122
0.0971817
0.9225821
0.9949371
Xkr4|Xkr4
4459.42848
0.6341912
0.2902594
2.1849118
0.0288953
0.2180911
Gm39585|Gm39585
5797.38480
1.1798031
0.2580070
4.5727554
0.0000048
0.0001973
Tcea1|Tcea1
2657.79689
0.1319203
0.3560044
0.3705581
0.7109667
0.9583981
Rgs20|Rgs20
422.79650
0.3482376
0.8144256
0.4275868
0.6689520
0.9469053
Gm16041|Gm16041
21.30106
7.1396504
3.9105846
1.8257246
0.0678918
0.3580246
Atp6v1h|Atp6v1h
10659.10495
0.1185235
0.3172387
0.3736100
0.7086945
0.9576220
LOC108167788|LOC108167788
11.33585
-0.5355535
3.7403736
-0.1431818
0.8861466
0.9949371
Npbwr1|Npbwr1
924.87540
-0.6449281
0.7097813
-0.9086292
0.3635459
0.7811879
Oprk1|Oprk1
15661.74234
0.2927281
0.3454541
0.8473720
0.3967878
0.8053124
St18|St18
279.69289
-1.0473877
0.9882176
-1.0598756
0.2892012
0.7120102
Rb1cc1|Rb1cc1
12217.86115
0.3066224
0.3041177
1.0082359
0.3133412
0.7348823
4732440D04Rik|4732440D04Rik
5369.16278
-0.1463996
0.2594995
-0.5641616
0.5726441
0.9132695
Gm26901|Gm26901
163.46252
-0.7535570
1.0236447
-0.7361509
0.4616389
0.8511309
Pcmtd1|Pcmtd1
4321.05905
0.5461501
0.2827538
1.9315391
0.0534164
0.3103856
Rrs1|Rrs1
625.86245
-0.1126443
0.7061219
-0.1595253
0.8732551
0.9949371
Mybl1|Mybl1
184.22714
-2.0915101
1.2989231
-1.6101878
0.1073569
0.4493503
Adhfe1|Adhfe1
452.68051
-1.0496649
0.8432838
-1.2447350
0.2132292
0.6256462
out of 23922 with nonzero total read count adjusted p-value < 0.1
LFC > 0 (up) : 628, 2.6% LFC < 0 (down) : 1046, 4.4% outliers [1]
: 2272, 9.5% low counts [2] : 2713, 11% (mean count < 6) [1] see
‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’
argument of ?results
1.1.3 Craving vs
Maintenance
Craving vs Maintenance CIVSA Results Overview
baseMean
log2FoldChange
lfcSE
stat
pvalue
padj
Mrpl15|Mrpl15
2162.26767
-0.2627261
0.3466449
-0.7579115
0.4485040
0.8651323
Gm39587|Gm39587
20.30129
-0.3516955
3.3976607
-0.1035111
0.9175574
0.9939251
Lypla1|Lypla1
994.90242
-0.0715356
0.8799222
-0.0812976
0.9352053
0.9951721
Xkr4|Xkr4
4459.42848
-0.6679818
0.2714805
-2.4605151
0.0138738
0.1495170
Gm39585|Gm39585
5797.38480
-1.0936780
0.2412810
-4.5327985
0.0000058
0.0002749
Tcea1|Tcea1
2657.79689
-0.5356827
0.3330620
-1.6083574
0.1077569
0.4613839
Rgs20|Rgs20
422.79650
-0.5371710
0.7618243
-0.7051114
0.4807409
0.8833096
Gm16041|Gm16041
21.30106
0.4549475
3.4906011
0.1303350
0.8963014
0.9909408
Atp6v1h|Atp6v1h
10659.10495
-0.6141458
0.2967693
-2.0694383
0.0385050
0.2733236
LOC108167788|LOC108167788
11.33585
-0.2865007
3.5020543
-0.0818093
0.9347983
0.9951721
Npbwr1|Npbwr1
924.87540
0.3369643
0.6640179
0.5074627
0.6118302
0.9372345
Oprk1|Oprk1
15661.74234
-0.4369729
0.3231413
-1.3522657
0.1762903
0.5932565
St18|St18
279.69289
-0.9608973
0.9259208
-1.0377748
0.2993749
0.7560424
Rb1cc1|Rb1cc1
12217.86115
-0.3863013
0.2844705
-1.3579665
0.1744743
0.5910084
4732440D04Rik|4732440D04Rik
5369.16278
-0.3454179
0.2427896
-1.4227046
0.1548218
0.5589467
Gm26901|Gm26901
163.46252
0.5242196
0.9577914
0.5473213
0.5841580
0.9274722
Pcmtd1|Pcmtd1
4321.05905
-0.3426769
0.2644209
-1.2959526
0.1949918
0.6240924
Rrs1|Rrs1
625.86245
-0.3724096
0.6606702
-0.5636845
0.5729689
0.9244236
Mybl1|Mybl1
184.22714
3.5366832
1.2149336
2.9110094
0.0036026
0.0603776
Adhfe1|Adhfe1
452.68051
0.8959769
0.7889411
1.1356702
0.2560946
0.7074992
out of 23922 with nonzero total read count adjusted p-value < 0.1
LFC > 0 (up) : 1084, 4.5% LFC < 0 (down) : 324, 1.4% outliers [1]
: 2271, 9.5% low counts [2] : 2713, 11% (mean count < 6) [1] see
‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’
argument of ?results
1.2
Annotate Results Files and Pull Significant DEGs
This subsection is for annotation the results files to include all
necessary information for downstream analyses.
After annotation I subset the results to only view the significantly
differentially expressed genes with an adjusted p-value > 0.05.
par-3 family cell polarity regulator beta, opposite strand 2
105243940
ENSMUSG0….
Col4a3|Col4a3
29.96706
19.601203
3.5629545
5.501390
0.0000000
0.0000079
Col4a3
collagen, type IV, alpha 3
12828
ENSMUSG0….
Dnajb3|Dnajb3
70.36681
9.093523
2.0877638
4.355628
0.0000133
0.0017489
Dnajb3
DnaJ heat shock protein family (Hsp40) member B3
15504
ENSMUSG0….
Asb18|Asb18
12.85015
18.822115
3.6835623
5.109759
0.0000003
0.0000497
Asb18
ankyrin repeat and SOCS box-containing 18
208372
ENSMUSG0….
Gm41920|Gm41920
91.39199
20.766926
3.6195797
5.737386
0.0000000
0.0000030
Gm41920
predicted gene, 41920
105246666
NA
LOC105246710|LOC105246710
41.29964
21.372039
3.6789651
5.809253
0.0000000
0.0000022
LOC105246710
NA
NA
LOC108167732|LOC108167732
59.78860
21.982685
3.2366754
6.791749
0.0000000
0.0000000
LOC108167732
NA
NA
Fam129a|Fam129a
58.26340
9.542099
2.6179113
3.644928
0.0002675
0.0256134
Fam129a
NA
NA
Fmo1|Fmo1
51.26689
21.615385
2.7780474
7.780783
0.0000000
0.0000000
Fmo1
flavin containing monooxygenase 1
14261
ENSMUSG0….
Gorab|Gorab
666.84329
2.774594
0.5327144
5.208407
0.0000002
0.0000307
Gorab
golgin, RAB6-interacting
98376
ENSMUSG0….
Dnah14|Dnah14
41.60299
21.334681
3.4180513
6.241767
0.0000000
0.0000003
Dnah14
dynein, axonemal, heavy chain 14
240960
ENSMUSG0….
Atf3|Atf3
433.56152
-2.808330
0.7372279
-3.809311
0.0001394
0.0145805
Atf3
activating transcription factor 3
11910
ENSMUSG0….
B230208H11Rik|B230208H11Rik
48.28040
21.761939
2.7106546
8.028297
0.0000000
0.0000000
B230208H11Rik
RIKEN cDNA B230208H11 gene
320273
ENSMUSG0….
LOC108167797|LOC108167797
82.55280
9.456510
2.0988898
4.505482
0.0000066
0.0009335
LOC108167797
NA
NA
Fam26f|Fam26f
45.72461
21.035286
2.7780345
7.572003
0.0000000
0.0000000
Fam26f
NA
NA
Mical1|Mical1
109.45019
10.585487
2.6621850
3.976240
0.0000700
0.0077657
Mical1
microtubule associated monooxygenase, calponin and LIM domain containing
1
171580
ENSMUSG0….
1.2.2 Maintenance vs
FoodTrained
Significant DEGs Maintenance vs FoodTrained
baseMean
log2FoldChange
lfcSE
stat
pvalue
padj
symbol
description
entrez
ensembl
Gm39585|Gm39585
5797.38480
1.179803
0.2580070
4.572755
0.0000048
0.0001973
Gm39585
predicted gene, 39585
105243853
NA
Snhg6|Snhg6
522.81858
1.260214
0.4234933
2.975758
0.0029227
0.0446701
Snhg6
small nucleolar RNA host gene 6
73824
ENSMUSG0….
Trpa1|Trpa1
47.61966
-22.462571
3.6001396
-6.239361
0.0000000
0.0000001
Trpa1
transient receptor potential cation channel, subfamily A, member 1
277328
ENSMUSG0….
Rpl7|Rpl7
2486.96588
1.005240
0.3038678
3.308150
0.0009391
0.0184408
Rpl7
ribosomal protein L7
19989
ENSMUSG0….
Gm28836|Gm28836
24.57062
20.693225
3.9082681
5.294730
0.0000001
0.0000064
Gm28836
predicted gene 28836
102635593
NA
Gm32790|Gm32790
44.07127
9.563380
2.2173604
4.312957
0.0000161
0.0005913
Gm32790
predicted gene, 32790
102635458
NA
Fer1l5|Fer1l5
230.33092
-2.687277
0.8236933
-3.262473
0.0011044
0.0209149
Fer1l5
fer-1-like 5 (C. elegans)
100534273
ENSMUSG0….
Sema4c|Sema4c
1298.41081
-1.822125
0.6102459
-2.985887
0.0028276
0.0434624
Sema4c
sema domain, immunoglobulin domain (Ig), transmembrane domain (TM) and
short cytoplasmic domain, (semaphorin) 4C
20353
ENSMUSG0….
Vwa3b|Vwa3b
68.09876
8.717777
2.3896076
3.648204
0.0002641
0.0065801
Vwa3b
von Willebrand factor A domain containing 3B
70853
ENSMUSG0….
4930439A04Rik|4930439A04Rik
239.71050
7.561383
1.1113671
6.803677
0.0000000
0.0000000
4930439A04Rik
RIKEN cDNA 4930439A04 gene
78119
NA
Il1r2|Il1r2
73.49166
9.935852
2.8549805
3.480182
0.0005011
0.0112029
Il1r2
interleukin 1 receptor, type II
16178
ENSMUSG0….
Il1rl2|Il1rl2
153.74171
-24.241136
3.7848523
-6.404777
0.0000000
0.0000000
Il1rl2
interleukin 1 receptor-like 2
107527
ENSMUSG0….
LOC108167622|LOC108167622
55.54376
21.760261
2.9333442
7.418244
0.0000000
0.0000000
LOC108167622
NA
NA
Gm10561|Gm10561
31.31746
-8.823060
2.8132289
-3.136275
0.0017111
0.0295377
Gm10561
predicted gene 10561
628004
NA
Clk1|Clk1
11902.12420
1.056339
0.3336618
3.165897
0.0015461
0.0272603
Clk1
CDC-like kinase 1
12747
ENSMUSG0….
Stradb|Stradb
1233.07664
1.338207
0.3402418
3.933107
0.0000839
0.0024506
Stradb
STE20-related kinase adaptor beta
227154
ENSMUSG0….
Fam117b|Fam117b
1077.36822
1.740967
0.5934892
2.933444
0.0033522
0.0491723
Fam117b
family with sequence similarity 117, member B
72750
ENSMUSG0….
D230017M19Rik|D230017M19Rik
353.30341
1.817985
0.5121487
3.549721
0.0003856
0.0090159
D230017M19Rik
RIKEN cDNA D230017M19 gene
320933
NA
Arpc2|Arpc2
3532.07869
-1.115469
0.2391510
-4.664286
0.0000031
0.0001315
Arpc2
actin related protein 2/3 complex, subunit 2
76709
ENSMUSG0….
Cnppd1|Cnppd1
1597.07545
-1.205071
0.3232300
-3.728214
0.0001928
0.0049685
Cnppd1
cyclin Pas1/PHO80 domain containing 1
69171
ENSMUSG0….
1.2.3 Craving vs
Maintenance
Significant Craving vs Maintenance
baseMean
log2FoldChange
lfcSE
stat
pvalue
padj
symbol
description
entrez
ensembl
Gm39585|Gm39585
5797.38480
-1.0936780
0.2412810
-4.532799
0.0000058
0.0002749
Gm39585
predicted gene, 39585
105243853
NA
Trpa1|Trpa1
47.61966
20.0711689
3.3932805
5.914975
0.0000000
0.0000004
Trpa1
transient receptor potential cation channel, subfamily A, member 1
277328
ENSMUSG0….
Rpl7|Rpl7
2486.96588
-1.0763933
0.2841797
-3.787721
0.0001520
0.0049557
Rpl7
ribosomal protein L7
19989
ENSMUSG0….
1700001G17Rik|1700001G17Rik
704.98525
1.5933623
0.4255143
3.744557
0.0001807
0.0056755
1700001G17Rik
RIKEN cDNA 1700001G17 gene
67503
ENSMUSG0….
Fer1l5|Fer1l5
230.33092
2.7365069
0.7714625
3.547168
0.0003894
0.0105650
Fer1l5
fer-1-like 5 (C. elegans)
100534273
ENSMUSG0….
Sema4c|Sema4c
1298.41081
1.9826591
0.5709143
3.472779
0.0005151
0.0131477
Sema4c
sema domain, immunoglobulin domain (Ig), transmembrane domain (TM) and
short cytoplasmic domain, (semaphorin) 4C
20353
ENSMUSG0….
Il1rl2|Il1rl2
153.74171
22.2736364
3.5649386
6.247972
0.0000000
0.0000001
Il1rl2
interleukin 1 receptor-like 2
107527
ENSMUSG0….
Fam126b|Fam126b
8348.97120
-0.8264248
0.2373035
-3.482565
0.0004966
0.0127789
Fam126b
NA
NA
Pard3bos2|Pard3bos2
48.12983
-23.7390363
3.6316177
-6.536766
0.0000000
0.0000000
Pard3bos2
par-3 family cell polarity regulator beta, opposite strand 2
105243940
ENSMUSG0….
Zdbf2|Zdbf2
3291.47772
-0.9550029
0.3014901
-3.167610
0.0015370
0.0313656
Zdbf2
zinc finger, DBF-type containing 2
73884
ENSMUSG0….
Klf7|Klf7
3114.99312
-0.6650793
0.2115799
-3.143396
0.0016700
0.0335025
Klf7
Kruppel-like factor 7 (ubiquitous)
93691
ENSMUSG0….
Map2|Map2
15113.61849
-0.6774675
0.1781094
-3.803658
0.0001426
0.0047204
Map2
microtubule-associated protein 2
17756
ENSMUSG0….
Gm31047|Gm31047
43.06422
19.1697461
3.6817645
5.206674
0.0000002
0.0000130
Gm31047
predicted gene, 31047
102633145
NA
Arpc2|Arpc2
3532.07869
0.7956339
0.2237915
3.555246
0.0003776
0.0103047
Arpc2
actin related protein 2/3 complex, subunit 2
76709
ENSMUSG0….
Cnppd1|Cnppd1
1597.07545
0.9103462
0.3024970
3.009439
0.0026173
0.0473299
Cnppd1
cyclin Pas1/PHO80 domain containing 1
69171
ENSMUSG0….
Dnajb2|Dnajb2
1487.62355
1.1910673
0.3095799
3.847366
0.0001194
0.0040449
Dnajb2
DnaJ heat shock protein family (Hsp40) member B2
56812
ENSMUSG0….
Atg9a|Atg9a
8813.08524
2.0179237
0.4649911
4.339704
0.0000143
0.0006226
Atg9a
autophagy related 9A
245860
ENSMUSG0….
Resp18|Resp18
56768.10866
-0.8522294
0.2302091
-3.701980
0.0002139
0.0065238
Resp18
regulated endocrine-specific protein 18
19711
ENSMUSG0….
Chpf|Chpf
11304.97381
2.0829961
0.3845635
5.416520
0.0000001
0.0000049
Chpf
chondroitin polymerizing factor
74241
ENSMUSG0….
Tmem198|Tmem198
3311.19558
1.5407766
0.3702558
4.161384
0.0000316
0.0012693
Tmem198
transmembrane protein 198
319998
ENSMUSG0….
1.3 Visualize Shrinkage
Estimations
Following subsetting, I want to see how the data looks based on
effect size so I have run a shrinkage model using the “apeglm” method
which normalizes counts based on effect size without major changes to
data structure. I have included non-shrunk models as well for
comparison.
1.3.1
Non-Shrunken L2FC
1.3.1.1 Craving vs
FoodTrained
Version
Author
Date
dc04b15
avm27
2024-07-17
1.3.1.2 Maintenance vs
FoodTrained
Version
Author
Date
dc04b15
avm27
2024-07-17
1.3.1.3 Craving vs
Maintenance
Version
Author
Date
dc04b15
avm27
2024-07-17
1.3.2
Shrunken L2FC
1.3.2.1 Craving vs
FoodTrained
Version
Author
Date
dc04b15
avm27
2024-07-17
1.3.2.2 Maintenance vs
FoodTrained
Version
Author
Date
dc04b15
avm27
2024-07-17
1.3.2.3 Craving vs
Maintenance
Version
Author
Date
dc04b15
avm27
2024-07-17
2 Looking at
Differentially Expressed Genes
We can generate a matrix to identify where our differentially
expressed genes lie between each group based on an alpha value of
0.05.
3.1
Normalized Counts and Plots of Genes of Interest
level_order <-c("FoodTrained", "Maintenance", "Craving") # this vector might be useful for other plots/analysesnormalized_counts <-counts(dds, normalized =TRUE)
Volcano plots are generated by gathering the FDR corrected p-values
from each analysis. The adjusted p-values undergo a -log10
transformation to generate these FDR corrected values. Any rows
containing “NA” are ommitted from analysis. Data points are colored
based on increasing or decreasing value and plotted using ggplot2.
Adjusted p-value cutoff is 1.3 following log transformation.(-log10 0.05
= 1.3)
Interactive Volcano plots are generated by Using the same method as
above, but are plotted so that genes may be identified and log2FC can be
seen from each gene. Plots can be used to customize which genes appear
based on what you would like to see.
4 Pathway Analysis
4.1
Pathway Analysis GO
First, the background for the GO universe must be set. In order to do
so we will subset the first results file “res” for any genes that map to
entrez IDs and then we only take the unique ones. This leaves us with
~31,000 genes for the remaining background which we list as our midbrain
dopamine neuronal expressed genes.
Pathway analysis is conducted using enrichGo and enrichKEGG based on
a p-value of 0.05 and a q-value of 0.10 and including all significantly
differentially expressed genes.
# Remove any genes that have greater than 1.5 LFC Standard Errorresults_sig_entrez <-subset(results_sig, lfcSE <=1.5)results_sig_entrez2 <-subset(results_sig2, lfcSE <=1.5)results_sig_entrez3 <-subset(results_sig3, lfcSE <=1.5)# Remove any genes that have greater than adj.pvalue 0.05results_sig_entrez <-subset(results_sig_entrez, padj <=0.05)results_sig_entrez2 <-subset(results_sig_entrez2, padj <=0.05)results_sig_entrez3 <-subset(results_sig_entrez3, padj <=0.05)# Remove any genes that do not have any entrez identifiersresults_sig_entrez <-subset(results_sig_entrez, is.na(entrez) ==FALSE)results_sig_entrez2 <-subset(results_sig_entrez2, is.na(entrez) ==FALSE)results_sig_entrez3 <-subset(results_sig_entrez3, is.na(entrez) ==FALSE)# Create a matrix of gene log2 fold changesgene_matrix <- results_sig_entrez$log2FoldChangegene_matrix2 <- results_sig_entrez2$log2FoldChangegene_matrix3 <- results_sig_entrez3$log2FoldChange# Add the entrezID's as names for each logFC entrynames(gene_matrix) <- results_sig_entrez$entreznames(gene_matrix2) <- results_sig_entrez2$entreznames(gene_matrix3) <- results_sig_entrez3$entrezgo <-enrichGO(gene =names(gene_matrix),OrgDb ="org.Mm.eg.db",readable = T,universe = gene_matrix_background,ont ="ALL",pAdjustMethod ="fdr",pvalueCutoff =0.05,qvalueCutoff =0.10)go2 <-enrichGO(gene =names(gene_matrix2),OrgDb ="org.Mm.eg.db",readable = T,universe = gene_matrix_background,ont ="ALL",pAdjustMethod ="fdr",pvalueCutoff =0.05,qvalueCutoff =0.10)go3 <-enrichGO(gene =names(gene_matrix3),OrgDb ="org.Mm.eg.db",readable = T,universe = gene_matrix_background,ont ="ALL",pAdjustMethod ="fdr",pvalueCutoff =0.05,qvalueCutoff =0.10)
4.1.1 Craving vs
FoodTrained
kable(head(go, 10), caption ="Top 10 GO Pathways for Craving vs FoodTrained", format ="html", align ="l", booktabs =TRUE) %>%kable_styling(html_font ="Arial", font_size =10, bootstrap_options =c("striped", "hover", "condensed", "responsive")) %>%scroll_box(width ="1000px", height ="500px")